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Abstract 

The problem consists of an electromagnetic wave incident on one of the faces of a rectangular 
waveguide. Part of the wave is coupled to the waveguide and part of it propagates outside, assumed 
to be free space. The shape of the propagating wavefront is to be found out, both inside and outside 
the waveguide. The electric field propagation can be found out by solving Maxwell's equations both 
inside and outside the waveguide, with proper boundary conditions. The numerical technique used 
here to solve Maxwell's equations is Finite Difference Time Domain method, which utilizes central 
difference approximation. Two different boundary conditions are imposed on the waveguide and 
on the free space boundary. Waveguide walls are characterised by PEC and free space boundary is 
characterised by PML boundary conditions. The code is written in 'C language, output is plotted 
in MATLAB.An mpeg movie is created from the code generated data, showing the evolution of 
the electric field with time inside the waveguide 



MATHEMATICAL FORMULATION 

The mathematical formulation of the problem consists of the specification of the waveg- 
uide parameters,imposing the proper boundary conditions, applying central difference ap- 
proximation to the Maxwell's equations and discreetizing them in all the 3 co-ordinates, which 
would be used in the program. 

waveguide specification 

The dimension of the waveguid is a=2.29 cm, b= 1.02cm. The frequency of the incident 
wave is 10 Ghz. The waveguide length is 5 x A,where A is the wavelength of the incident 
wave which in this case is 3 cm. 



Boundary Conditions 

The problem spans in two mediums, inside the waveguide and outside i.e in free space. The 
total space including these regions are discreetized into grids or cells,with different boundary 
conditions for both the regions. The number of cells per wavelength is 25. The boundary 
surface of the waveguide is a perfect electrical conductor. PEC boundary condition is 

1. Tangential component of the electric field is in the waveguide at the boundary. 

2. En is not necessarily zero in the wave guide at the boundary as there may be surface 
charges on the conducting walls 

The continuity of the magnetic field implies that normal component of magnetic field is 
at the boundary. There are currents induced in the guides but for perfect conductors these 
can be only surface currents. Hence, there is no continuity for H t . 



Far field conditions 

The far field conditions must be satisfied at the end of the vaccum grid. Here perfectly 
matched layers are imposed as boundary conditions. These represent the absorbing boundary 




FIG. 1. one single fdtd cell geometry 



condition that the field must vanish at infinity. But numericaly no grid can extend upto 

infinity.The last few surfaces of the vaccum cell are imposed with PML conditions. A PML 

is designed to absorb waves coming to its direction,Wherever an x derivative d/dx appears 

in the wave equation, it is replaced by: 

d 1 d 

dx 
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Where u) is theangular frequency and a is some function of x. Wherever a is positive, propagating 
waves are attenuated, because 



k 
exp[i(kx — u)t)] — > exp[i(kx — cot) / a(x')dx' 



(2) 



This transformation of coordinates causes waves to attenuate whenever they have a depen- 
dance of expiikx). 

source implementation 

The length of the guide is along x direction, the width along y direction and the height 
along z direction. Plane wave emanates from the x = face. Thus all the E z components on 
the x = face are taken as the stimulus. The source is taken as a gaussian pulse expressed 

as 

E z (i,j, k)(t) = Aexp[-c(t - i /O 2 ] (3) 

Where c is a constant. 

Specification of cells 

The number of cells on each direction i.ex,y,z are taken as the same,here it is 50. The first 
25 cells on each direction corressponds to the inside of the waveguide and the outermost 
of these cells on each direction provide the PEC boundary. The next 20 cells represent the 
free space grid and the last 5 cells on each direction represents the perfectly matched layers, 
where the wave gradually decays without any reflection. Say, rti is the total number of cells 
in each direction. (z=a;,t/, z), then the space differntial i.e di is given by 

dy = guidewidth/riy (4) 

dz = guide height /n z (5) 

dx is chosen to be the smallest of either dy or dz. The time step is defined as 

dt = l/(c x sqrt[(l/dx 2 ) + (1/dy 2 ) + (1/dz 2 )]) (6) 

Total stimulated time is given by 

T = maxiteration x dt (7) 
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memory allocation 

Mesh is set up such that the tangential E vectors form the outer faces of the simulation 
volume. Say there are n x ,n y , n z cells in the mesh, but the number of E x vector components 
are given by 

n x (n y + l)(n z + 1) (8) 

The number of E y components are given by 

(n x + l)n y (n z + 1) (9) 

similarly for E z components also. 

The H arrays are staggerred half a step off from the E arrays in every direction, thus H 

arrays are one cell smaller than E cells in each direction. 

Thus outer faces of the mesh consist of E componnts only and the inner mesh consists of 

the H components. 



Update equations 

For a source free medium (vaccum) the maxwell's equations are 

f^V** (10) 

dt ix Q 

Where E and H are vectors in 3 Dimensions. These equations interms of the x,y,z compo- 
nents can be written as 

dE x l r dH z 8H V ^ 

~9f = -[-*r - -*r] ( 12 ) 

dEy 
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dE z 

dt eo L dx dy i 
The equations for the H vector interms of its components are 
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(13) 
(14) 

(15) 



dH z = 1 dE y dE 

dt /i dx dy 

Equations5tol0can be interpreted as a small change in E x over time is related to a 

change mH y SindH z m space. Finite Difference approximation can be used for the above 

equations. The resulting equations can be expressed as 

E x (z',t> + f)- E x (z',t> -f)l My' + f ,t') - H z {y> - f ,t') 



At e At/ 

H y (z' + f,t')-H y (z'-^,t') . 
Az 

E y (z',t' + %) - E y (z',t' - f ) 1 H x (z' + A£,t') - H x {z' - *f,t) 



(18) 



At e Az 

H z (x' + ^,t')-H z (x'-^,t' 



Ax 
E z (z', t' + f ) - £ z (*', f - %) 1 #,(x' + ^, t') - /^(x' - A£, t') 



(19) 



At e Ax 



(20) 
Ay 

Similar equations can be written for H field with central difference approximation. Alittle 
manipulation gives the update equations used in the code. These are 

At At At H (z' + — t') - H (z' - — t') 

Ex ( z ',t'+— ) = E x (z',t' -—) + —[ x[ 2 ' ' 2 ' j 



2 ' v ' 2 ' e A2 

#,(*' + ^,f)-g,(s'-^,0 

Ax 

F ,, fl | At _ At At /j^ + f ,f) - g^ - f ,f) 

tf 2 (x' + f:,t')-/^(x'-AE,t'). 



(21) 



Ax 



E z (z',t> + ^) = E z{ z',t> ^ ■ At^ + f,t')-^-f,0 



(22) 



1 ' 2 7 e Ax 

H x (y'+^,t')-H x (y'-^,t' 



(23) 
Ay V ; 

These equations are directly used to generate volumetric data for the simulation of the 
electromagnetic wave. 



LOGIC OF THE PROGRAM 

The logic is defined as: 

1. Define the parameters of the model waveguide. Discreetize the total space into grids. Specify 
the number of cells in each direction. Define the space step. 

2. Define the time step. Allocate dynamic memory for the field values to be calculated in 
each iteration. 

3. Set the boundary conditions and the source. 

4. For each time step, using the update equations for both the E and H fields, the prop- 
agation of the electromagnetic field is calculated through the whole grid, this process 
is repeated for each time step. 

5. All the data for this time and space stepping is outputed in .txt files and each file, when 
plotted,gives a picture of the E field propagation. When all these files are poltted 
sequentialya picture of the wavefront propagating through the waveguide is visualised. 

Complexity of the program 

The program uses the update equations as the basic equations for calculations. Equations 
21 to 23 are the update equations. During each of the time step,the field values are calculated 
n x xn y xn z , n x ,n y ,n z are the number of cells in the x, y, z directions respectively. The output 
of the program gives the E z value only.For this, the complexity of the algorithm is of the 
order n 3 where it is assumed that n x =n y =n z .The complexity for all time steps is n 3 x T 
where T is the total number of timesteps. 

Error propagation 

The error in this method comes from the fact that the total no of cells multiplied by the 
space steps must equal to the dimensions of the waveguide. In this case n y x dy = Guide 
width in meters, and n z x dz = guide height in meters, dx is chosen to be the smallest of 
dy,dz. Another source of error is the size of the time step. This time step size should be such 



that the plane wave travels a single cell length in a single time step. Otherwise FDTD can not 
keep up with signal propagation, because FDTD computes a cell only from its immediate 
neighbours. 

DATA AND ANALYSIS OF THE ACTUAL PROGRAM PERFORMANCE 

The program output was interms of 200 .txt files. These files contain 3D data for the E z 
field vector in each timestep. The performance of the program was listed in another .txt 
files whose contents are displayed below. 

program performance 

Dynamically allocated bytes : 1173832 bytes. 

Data and figures 

The program uses one thousand time steps. In each time step,the program calculates the 
E z values throughout the grid. Two hundred data files are generated to simulate the wave 
propagation through the waveguide. These files are plotted using MATLAB environment. An 
mpeg movie has also been created using the data files. 

Fig(2) 

Fig(2) is a plot of the simulated time against the number of space steps. Number of 
time steps is plotted along x axis and the simulation time is plotted along the y axis. The 
relationship is linear as shown in the graph. As the wave propagates inside the waveguide,the 
complexity increases and the system time also increases. Fig(2) supports this observation. 

Fig(3)and Fig(4) 

Fig(3) is a three dimensional plot of E z at a particular time step. The plot is generated 
using the data file for the 995 th time step. The software used for plotting was MATLAB. 
Fig (4) is a three dimensional plot of the propagating wavefront. 
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TABLE I. 

Below are the data for system time taken to compute the E z values in the grid in each time 
step. The space steps calculated during the program are dx = 0.00113333 meters, dy = 0.001145 

meter,dz = 0.00113333. The total simulation region is 0.0229 x 0.0102 x 0.149896 meter 3 . The 
total time for simulation is 2.19002 x 10 -09 seconds. The complete data is not mentioned here. 

time step system time 



1 


2.19002e-12sec 


2 


4.38005e-12sec 


3 


6.57007e-12sec 


4 


8.7601e-12sec 


5 


1.09501e-llsec 


6 


1.31401e-llsec 


7 


1.53302e-llsec 


8 


1.75202e-llsec 


9 


1.97102e-llsec 


422 


9.2419e-10 


423 


9.2638e-10 


424 


9.2857e-10 


425 


9.3076e-10 


426 


9.3295e-10 


427 


9.3514e-10 


428 


9.3733e-10 


591 


1.2943e-09 


592 


1.29649e-09 


593 


1.29868e-09 


594 


1.30087e-09 


595 


1.30306e-09 


995 


2.17907e-09 


996 


2.18126e-09 


997 


2.18345e-09 


998 


9 2.18564e-09 


999 


2.18783e-09 
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FIG. 2. simulation time vs time step number 



A hollow rectangular waveguide can propagate TM and TE modes but not TEM 
modes. The TE modes are characterised by fields with E z =0 while H z must satisfy the 
wave equation 

o l x o l x 
Similarly the E z fields in the TM modes satisfy the wave equation 

d 2 d 2 



o l x o z x 



(25) 



Where 



E z (x,y,z) = e z (x,y)exp{-jf3z) 



(26) 



and 



H z (x,y,z) = h z (x,y)exp(-j(3z) 



(27) 
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FIG. 3. shape of the wavefront in one time step 

k 2 c = k 2 -p 2 (28) 

k c is the cut off wave number. The general solution can be found out for E z by using separation 
of variables. We use 

e z (x,y)=X(x)Y(y) (29) 

The general solution for this equation is 



e z {x, y) = (Acosk x x + Bsink x x)(Ccosk y y + Dsink y y) 

The boundary conditions are directly applied to the e z component which are 

e z (x,y) = 

at x = 0, a. 

e z (x,y) = 
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(30) 



(31) 



(32) 



fdtdSOatxt" u 1:2:3 
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FIG. 4. Propagation of the wavefront 

aty = 0, b. 

Applying 31 and 32 we can get the solution that A = and k = rmr/a for m = 0, 1, 2, 3.... 
Also we get 

C = and k y = nir/b for n = 0, 1, 2, 3.... 
The solution for £^ can thus be written as 



mux niry 

E z {x,y,z) = B m nsin sin^ — exp{—jpz) 

a b 



(33) 



The transverse field components for the TM m n mode can be written as 



2 mirx niry 

E x = jrac—jprrmakB m ncos sin— — expi—jpz) 

a b 



(34) 
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2 mirx niry 
E y = frac—j prmbk c B m nsin cos——exp(—jpz) (35) 



H x = fracjuemrbklB m nsin cos — — exp(—jj3z) (36) 



—jcuermr rmrx wny 

H v = — — B m ncos sin— — expi—jpz) (37) 

akt a b 



From equations 34 to 36 we can say that these expressions are identicaly zero if either of m 
or n is zero. So the lowest order TM mode possible is TM\\. Its cut off frequency is 

W1 =2^r«< )2 +<I )2)) < 38) 

For this problem a = 2.29cm and b = 1.02cm. Thus the cutoff frequency for TM\\ mode is 
16.09 GHz.TMil mode is calculated because the program outputs only E z values which are 
possible for TM\\ modes. 

TEST FOR THE CORRECTNESS OF THE SOLUTION 

Figure 5 shows the electric field distribution in the waveguide for the TM\\ mode.Fig(4) 
showed a similar nature of propagation of the electric field vectors. From that it can be said 
that the program was able simulate the propagation of one of the eigenmodes, i.e the TM\\ 
mode. 

Although the free space was not simulated in this program,so the propagation of the wave 
outside the guide was not found out, so there is further room for improvement in this 
program. 

REFERENCE 

1. The waveguide specificationwas taken from David K. Cheng, Field and Wave Electro- 
magnetics, 2nd ed., pages 554-555. 
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FIG. 5. E field distribution in a rectangular waveguide 

2. A part of the program was taken from the website www.toyfdtd.com 

3. Part of the algorithm and workings of FDTD was taken from the webpage of EPSRC 
Summer School 2007, Ian Drumm. 

4. The mathematical analysis of the eigenmodes of propagation inside a waveguide was 
developed with help of the book Microwave engineering, David M.Pozar, Wiley publi- 
cations, Third edition. 
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